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1 INTRODUCTION 
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Stellar disks are one of the most conspicuous features of 
galaxies: they are the sites of quiescent star formation and 
display such distinctive features as spiral arms and bars. 
Current data suggest that 40-60% of all luminosity in the 
local Universe is contributed by stars located in galaxy disks 
(Tasca & White 2005; Benson et al. 2007). In cold dark mat- 
ter universes, where structure forms hierarchically, it is be- 
lieved that galactic disks are the first components to form, 
with elliptical galaxies and bulges of spiral galaxies form- 
ing later through the merging of pre-existing galaxy disks 
(Kauffmann, White & Guiderdoni 1993; Baugh et al. 1996). 
Understanding in detail how galaxy disks form and evolve is 
therefore crucial to our understanding of galaxy formation 
as a whole. 

Traditionally, semi-analytic models of galaxy formation 
have estimated the total masses of stars, gas, and metals con- 
tained in each galaxy disk, but not the radial distribution 
of these components (see Kauffmann et al. 1999; Somerville 
& Primack 1999; Cole et al. 2000; Hatton et al. 2003 for 
recent examples). Instead, they explicitly assume that the 
surface density has an exponentially declining radial depen- 
dence. Consequently, such models are forced to adopt highly 
simplified prescriptions for star formation, which are unable 
to capture the physical mechanisms involved in this complex 
process. 

Conversely, models of individual disks (e.g. Mo, Mao & 
White 1998a; Mo, Mao & White 1998b; Firmani & Avila- 
Reese 2000; Efstathiou 2000; van den Bosch 2001; Monaco 



ABSTRACT 

The formation of galactic disks and the efficiency of star formation within them are 
issues central to our understanding of galaxy formation. We have developed a detailed 
and versatile model of disk formation which combines the strengths of previous studies 
of isolated disks with those of hierarchical galaxy formation models. Disk structure is 
inferred from the distribution of angular momentum in hot halo gas and the hierar- 
chical build-up of dark matter, leading to theoretically generated systems where the 
evolution of surface density, rotation, velocity dispersion, stability and metallicity is 
predicted for annular regions of width 20-100pc. The model will be used to estab- 
lish whether the accepted theory of large scale structure formation in the Universe is 
consistent with observed trends in the properties of disk galaxies. 

This first paper explicitly examines the importance of embedding such calculations 
within a merging hierarchy of dark matter halos, finding that this leads to dramatically 
different formation histories compared to models in which disks grow in isolation. 
Different models of star formation are explored, and are found to have only a secondary 
influence on the properties of the resulting galaxy disks, the main governing factor 
being the inf ailing gas supply from the hot halo. 

2004; Kampakoglou & Silk 2006) study the radial distribu- 
tion of these properties in great detail but do not always 
consider the dynamical evolution of the system, digressions 
from the steady state or, perhaps more importantly, inter- 
actions with neighbouring galaxies. 



A similar dichotomy has existed within simulations of 
galaxy formation, which have often focused on either the de- 
tailed evolution of a single halo (Abadi et al. 2003) or on the 
properties of large samples of galaxies using comparatively 
low resolution (Navarro &: Steinmetz 2000). This situation 
can be addressed by introducing methods to increase reso- 
lution in areas of interest, beginning on cosmological scales 
(~ lOOMpc) but selecting smaller regions (~ 5kpc) to study 
using a softening lengths as small as 500pc (Sommer-Larsen, 
Gotz & Portinari 2003; Governato et al. 2004; Robertson et 
al. 2004). 

In order to make progress, the respective strengths of all 
these approaches must be brought together. Given sufficient 
information about the halos in any analytic or numerical 
model of large-scale structure, the profile of the disk which 
forms within each sub-halo can be determined with simple 
assumptions for the symmetry of the system and the conser- 
vation (or otherwise) of the angular momentum of halo gas 
as it cools and condenses to form a disk. Local values of the 
star formation rate, the rate of outflow of gas, metallicity 
and the velocity dispersion of gas clouds and so forth can 
then be estimated as a function of radius in the disk. 

This approach immediately provides three distinct ad- 
vantages over previous work. Firstly, it allows the surface 
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density profiles of galaxy disks to be determined from more 
fundamental premises (the distribution of halo gas, angular 
momentum conservation, the frequency of galaxy mergers 
and so on) rather than simply by assumption. The metalli- 
cies of stellar and gaseous components can be followed sep- 
arately, rather than assuming that they trace the overall 
density. Secondly, knowledge of the surface density profile 
and associated rotation curve and stability allows for the 
inclusion of star formation models based upon local density, 
local stability on so on. Finally, all of this can be modelled 
within a realistic hierarchical galaxy formation scenario. 

We have developed such an improvement within the 
Galform semi-analytic model (Cole et al. 2000) and the 
remainder of this paper describes the techniques employed 
(Sj2| and presents basic results (Sj3|.The effects of hierarchical 
structure formation and galaxy mergers on disk formation 
are studied using three alternative models for star formation. 
These comparisons are carried out initially in the simplified 
context of an isolated, individual dark matter halo and then 
in the context of a similar system formed within a fully hier- 
archical model of galaxy formation. Conclusions are drawn 

in g 



2 THE MODEL 

This section describes the methods used to compute the ra- 
dial structure and star formation within galaxy disks. The 
model can be embedded within the Galform semi-analytic 
model of galaxy formation — which is not described in detail 
here — although it can also be run independently of Gal- 
form. The reader is referred to Cole et al. (2000) and 
Baugh et al. (2005) for complete details of Galform. Only 
those parts of the model directly affected by recent modifi- 
cations are detailed here. Throughout this work, a standard 
ACDM model of cosmology is assumed, with values chosen 
for consistency with the latter publication (S2m = 0.3, J1a = 
0.7, fib = 0.04, a 8 = 0.93, ho = 0.7). 



model this compression by assuming that dark matter par- 
ticles conserve the adiabatic invariant quantity j during the 
process of baryonic collapse and that there is no shell cross- 
ing (dM(j)/dt = 0). Throughout this paper is that M(j) 
and M(r) always refer to the mass of material with radius 
less than r or specific angular momentum less than j. 

This calculation uses j = %/ GMr for the initial dark 
matter distribution. Simulations of dark matter contraction 
on galactic scales (Gnedin et al. 2004) have indicated that 
this approximation overestimates the central density by not 
allowing for orbital eccentricity in the dark matter. More 
recently however, Dutton et. al. (2007) made additional cor- 
rections for the flatness of the disk and found that the com- 
bined effects lead to little net change from the initial, basic 
assumption. We therefore retain this simple approach. 

The baryonic component is assumed to be shock heated 
to the halo virial temperature during collapse. This hot gas 
is supported by thermal pressure and so, as it subsequently 
begins to cool, it can fall towards the halo centre. During this 
infall we assume that the collective drift angular momentum 
of the gas about some common axis is conserved (this angu- 
lar momentum arises as a result of tidal interactions between 
halos; Hoyle 1949). This claim has not always been sup- 
ported by simulations (Navarro & White 1994; Navarro & 
Steinmetz 2000) though recent work (Governato et al. 2004) 
indicates that this may well have been an issue of resolution. 
The assumption of conservation of angular momentum dur- 
ing collapse is one that can be changed within our model, a 
flexibility that will be used to advantage in future work. 

The distribution of density and angular momentum for 
simulated hot gas halos was recently found to closely match 
the following forms (Sharma & Steinmetz 2005): 



1 + 



dM(j) 

oc j exp 



or 



(2) 



2.1 Halo Structure 

In the ACDM cosmogony, the mass density of the Universe 
is dominated by cold dark matter, which interacts only via 
gravity (Bertone, Hooper & Silk 2005). Simulations of the 
behaviour of such particles (Navarro, Frenk, & White 1997) 
show that they undergo gravitational collapse into haloes 
with radial density dependance: 



The virial radius, r v , encloses the correct mean density for 
spherical top-hat collapse in the particular cosmology, thus 
specifying p v . For the cosmology used, this is 330 times the 
mean density at z = (Eke, Cole & Frenk 1996). The only 
free parameter, c, is found to be strongly correlated (and 
decreasing with) total halo mass (Navarro, Frenk, & White 
1995; Bullock et al. 2001), albeit with a large scatter. 

The pressureless dark matter particles are supported 
against further collapse only by their orbital motion. As 
baryonic material cools and condenses towards the centre 
of the halo, the dark matter density profile will be adiabat- 
ically compressed. Following Blumenthal et al. (1986), we 



Both expressions are normalized so that the appropriate to- 
tal baryonic mass, (fib/^M)Af v (where M v is the virial mass 
of the halo), is contained within the virial radius, as is the 
total angular momentum: 



Jtot = MvO') = AGA/ 5/2 /|£| 1/2 . (3) 



This defines the spin parameter, X. The value of this parame- 
ter, of a and th are particular to each halo, being distributed 
around typical values of 0.03, 0.9 and r v /20 respectively. 
Figure [I] shows the distribution of a found by Sharma & 
Steinmetz (2005). 

To relate the angular momentum distribution to the 
radially dependent density, we follow Efstathiou (2000) by 
assuming that cylindrical shells in the halo share the same 
j. Converseley, Firmani & Avila- Reese (2000) and van den 
Bosch (2001) both assume that the two components have 
the same density distribution at virialisation and that each 
spherical shell is in solid body rotation. 
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2.2 Disk Formation 

Though the details of non-linear gravitational collapse are 
most accurately modelled by N-body simulation, continu- 
ing to track the three dimensional location of the gas parti- 
cles as they cool is computationally impractical as part of a 
cosmological-scale model (the resolution achieved by Knebe, 
Dominguez & Dominguez-Tenreiro (2006), for example, is of 
the order of a kpc - too large to resolve individual galaxies 
accurately). By using the fact that the system has non-zero 
angular momentum and hence a well denned axis, the cool- 
ing phase can be modelled assuming cylindrical symmetry. 
This greatly reduces the complexity of further calculations 
without excessive departure from a realistic physical picture. 

We use the cooling model described in detail by Cole 
et al. (2000), which uses the cooling function A(T) = e/n 2 
calculated by Sutherland & Dopita (1993). Here T is the gas 
temperature, e is the energy loss rate per unit volume and 
n is the number density of hydrogen. The time taken for gas 
particles to cool is therefore 



icoolM = 



3 /imH ksT 
2 PgW A(T, zy 



(4) 



where gas density is p g , mean molecular mass is fimn and 
fee is Boltzmann's constant. As it cools, the gas will lose 
pressure, falling out of hydrostatic equilibrium and moving 
towards the centre of the halo until its residual angular mo- 
mentum lends sufficient rotational support against further 
collapse (Fall & Efstathiou 1980). The time taken for this to 
happen is assumed to be the greater of t coo i and the freefall 
time, 



ta(r) 



2GM(rO dr , 



-1/2 



dr ' 



(5) 



In the idealized situation that the angular momenta of all the 
gas shells are mutually aligned, the result will be a thin disk 
of gas with the original angular momentum vector normal 
to the disk plane. 

Both gas and dark matter components are assumed to 
conserve their specific angular momenta during the collapse 
process. Therefore, all baryonic material and dark matter 
which has the same initial specific axial angular momentum 
will share a common final radius in the plane which can be 
estimated b}{^] 

f 

R{j) = G[M diBk (j) + M DM 0')]' (6) 
where j is the specific angular momentum of circular orbits 
at radius R in the disk plane and Mi(j) is the mass of each 
component with specific angular momentum less than j. 

The disk is broken up into 200 concentric annuli be- 
tween the origin and the final radius of the most recently 
cooled gas. These are spaced linearly in specific angular mo- 
mentum. After each timestep, every annular region acquires 
all the halo gas within two cylindrical bounds, given by the 
range in specific angular momentum of that region, and two 
spherical bounds enclosing the portion of the halo which has 



1 An important underlying assumption is that the specific angu- 
lar momentum of the halo gas increases monotonically with its 
distance from the axis. 
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Figure 1. Values of the parameter a from a sample of 41 halos in 
a smoothed particle hydrodynamics (SPH) simulation (Sharma & 
Steinmetz 2005). This parameter controls the shape of the angu- 
lar momentum distribution throughout the halo, as described by 
eqn. pi. 



just cooled. After all masses have been updated, the function 
R(j) is reevaluated resulting in a revised profiles Mdisk(-R) 
and Mom(R)- The masses, mean metallicities and velocity 
dispersions of stars and gas can then be tracked for each of 
these regions. As such, individual annuli have widths of 20- 
lOOpc in typical disks. This choice is made so that regions 
are large enough for the use of scaling laws and locally av- 
eraged values to be appropriate, yet small enough so that a 
detailed picture of the disk profile is produced. Appendix |A"| 
shows that the final results are not particularly sensitive to 
this choice. 

Following thousands of disk profiles in such detail is 
computationally impractical and, as it turns out, largely un- 
necessary. The properties of each disk as a function of spe- 
cific angular momentum are therefore mapped onto a new 
grid of only 20 annuli after each calculation of gas infall 
and star formation. This coarser grid is used for the rest of 
the calculations (e.g. for computing galaxy properties dur- 



ing mergers, see £ 3.3 1 . Increasing this number above 20 pro- 



duces negligible change in the shape of the eventual profiles, 
as illustrated in Appendix [X] The full set of 200 annuli is 
then recreated by interpolation within this coarser grid to 
calculate gas cooling, star formation and outflow in each 
particular disk. 

Clearly eqn. (|6| incorrectly approximates the potential 
due to the disk as being spherically symmetric. The alterna- 
tive option, treating the disk as perfectly thin, is also only 
an approximation to reality and considerably complicates 
the calculations, mainly because there is no analytic form 
R(j,M). The relation between the three variables in this 
case would be (Cuddeford 1993): 



f{R) 



AGR 
2tt 



ada 



j,2 da 



dM(R') 



(7) 



An initial estimate of the mass profile M(R) then has to 
be corrected iteratively until this equation and Q are both 
satisfied. 

As shown in Fig. |2j there is surprisingly little differ- 
ence in results between the two approximations. Though 
they would differ by as much as 20% for an exponential 
disk, the overriding contribution of the dark matter and the 
relative central concentration of these disks reduce the typ- 
ical discrepancy to less than 5%. The less computationally 
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2 4 6 8 10 12 14 4 8 12 

r / kpc r / kpc 

Figure 2. Galactic disks which would form if all the gas in a 5 X 10 Mq virialised halo were to cool onto a plane, conserving the 
net component of its initial angular momentum as it does so. This constitutes a total baryonic mass of 7 X 10 10 Mq. in this example, 
the spin parameter, A = 0.03 and dark matter contraction parameter c = 10. The top left panel shows the mass of all the different 
components enclosed at each radius with surface density shown in the lower left panel. Circular velocity is plotted on the right. Original 
density and angular momentum profiles are given in equations JlJ and (J2J . Each line corresponds to a slightly different value of a (which 
governs the distribution of angular momentum in the original halo). The solid lines in the rotation curves plots were calculated using 
the approximation j 2 = GMR, while the dotted lines use equation tfl which approximates the disk as being perfectly thin. 



demanding approximation, eqn. ^6J, will therefore be used 
throughout this work. 

Figure [2] also shows the corresponding mass and surface 
density profiles of the disk (along with mass profiles for the 
original gas halo and dark matter) for each value of a. One 
clear feature to note at this stage is the sharp rise in surface 
density near the axis, leading to rotation curves with central 
peaks far greater than those observed, particularly for lower 
values of a. 

This feature is also seen in the model used by Dalcan- 
ton, Spergel & Summers (1997), which treats the halo as 
a uniformly rotating sphere. Given the simplicity of that 
assumption, the authors claimed no profound importance 
for these bulge-like centers. In the present case, the initial 
conditions are derived instead from simulations of halo for- 
mation and their lack of low-radius resolution provides a 
similar source of uncertainty. 

This problem of excess low angular momentum material 
persists even in simulations of disk systems (Abadi et al. 
2003) though, encouragingly, increased resolution does seem 



to lead to more observationally consistent mass distributions 
(Governato et al. 2006). 

Firmani & Avila-Reese (2000) and van den Bosch 
(2002) investigate this issue in detail, finding that systems 
with a lower baryon mass fraction (0.04-0.09) or higher spin 
(0.03-0.08) have more observationally consistent rotation 
curves. The factors responsible for the excessively peaked 
rotation curves here are therefore (1) a halo spin parameter 
of 0.03, (2) complete cooling of the baryonic component and 
(3) a universal baryon mass fraction of 0.13. 

Each of these differences are justifiable: Halos with 
A > 0.03 are rare in the simulations, so that value is appro- 
priate, van den Bosch (2002) finds that the fraction of the 
cooled baryonic mass approaches one for halo masses below 
~ 10 11 Mq/i _1 suggesting that complete cooling may not, in 
fact, be entirely unphysical for this halo. Lastly, the choice 
of cosmological parameters here actually underestimates the 
baryon fraction by comparison with the latest values from 
WMAP of fi M = 0.29 and fl b = 0.047 (Spergel et al. 2003). 
In conclusion, these choices are sufficiently well motivated 
to suggest that there may indeed be some missing physics. 
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One obvious candidate is gas outflow due to supernova 
winds but, when this is taken into account in ^3] this model 
suggests that it is not a significant enough effect to fully 
resolve the discrepancy. 

The assumption of a fully-formed initial halo will also 
exacerbate this effect. More realistic mass aggregation histo- 
ries are investigated by Firmani & Avila-Reese (2000) to the 
conclusion that, in earlier formed halos, gas cools at higher 
redshift when it is denser and consequently reaches equilib- 
rium at smaller radii. In the light of this, future discussions 
on dynamics will be restricted to the correct hierarchical 



case ([ 3.2 & { 3.3 1, leaving this isolated situation as merely 
a stage for introducing aspects of the disk's internal evolu- 
tion. 



2.3 Disk Evolution 

Once the disk has begun to form, the surface densities 
of gas, Eg and stars, E* in a given annular region are 
evolved using a set of coupled differential equations. These 
describe the star formation and feedback process whereby 
supernovae both expel gas from the surrounding region and 
heat any which remains. The former effect depletes the fuel 
supply from which stars are created while the latter can 
promote the efficiency with which they form (McKee & 
Ostriker 1977). This is therefore a complex loop, controlled 
ultimately by the velocity dispersion of the existing clouds 
within the disk, cr g , and the infall of new gas, Ei n (due 



to the cooling of halo gas described in £ 2.2 ) . The entire 
system of equations describing this system is detailed below. 
Where necessary, we give alternate forms for three different 
models of star formation and feedback. The models used are: 

(a) : The same fraction of gas is converted into stars every 
epicyclic period, a relationship which has empirical support 
(Kennicutt 1998); 

(b) : Star formation is promoted by gravitational instability, 
as proposed by Wang & Silk (1994); 

(c) : A model in which the star formation rate is propor- 
tional to the rate of gas cloud collisions. 



Star Formation Rate 



E* = e*— — / r msNiisN 



T* = < 





(8) 


o 


(a) 


Q 

Vi-Q 2 


(b) 


E E°"g 


(c) 



The parameter e* = 0.02 is the fraction of gas converted into 
stars per orbital period, r , as has been empirically deter- 
mined (Kennicutt 1998). Some fraction of this initial mass 
will eventually be returned to the ISM through supernova 
explosions and stellar winds. A Kennicutt IMF gives this 
fraction to be f x = 0.41 (Cole et al. 2000), where the mass 
of stars formed per supernova msN = 125M0. 

In model (b), the star formation rate is assumed to de- 
pend upon the epicyclic period, 



R 3 / 



dR 



1 1/2 



(9) 



and promoted by gravitational instability, as measured by 
the parameter 



Q 



TtGt c 



+ 



(10) 



This implies that there is no star formation in stable disks 
(Q > 1). The stellar velocity dispersion <r* = cr g /5. 

In model (c), the timescale for star formation, r*, is 
equal to the approximate time between collisions of molec- 
ular gas clouds, a possibility explored by Komugi, Sofue & 
Egusa (2007). Basic kinetic theory for clouds of mass m c i 
and radius r c i suggests 

•^mfp 



m c i 
r c iE E a E 



(11) 



The gas surface density is E g and A m f p is the mean free 
path of the gas clouds. The size and density of molecular 
gas clouds were found by Dame et al. (1986) to be related 
by pd ~ r ^i 1,3 which implies that the overall dependence 
on cloud properties in equation (111 will be relatively weak: 
E* oc r~®' 7 . With this in mind, and assuming that cloud sizes 



are fairly independent of the host system, eqn. (Ill leads to 
eqn. |8j;) by using characteristic values. The precise choice 
is not that important as a certain amount of uncertainty is 
implicitly contained in the free parameter e*. This work uses 
the value 



(Wei) = 2 x I^Mqpc- 1 



2 x 10 10 



years 



M pc 



In this theory, star formation would be promoted by a high 



velocity dispersion. Equation ( 10 1 can be rearranged to 
show that, in (b), it would be promoted by low velocity 
dispersion so the two models are therefore in direct contrast 
in this respect. The comparative effect this has on the 
evolution of the system is discussed in §3.1| 



Supernova Rate 



(12) 



C(r)E*(t-r)dr 



This calculates the current supernova rate per unit area, 
hsN, from the star formation history. The function £(r) is 
the initial mass function expressed in terms of the stellar 
lifetimes (rather then, as usual, their mass). The limit on 
the integral, r max ~ 50Myr, is therefore the lifetime of the 
lowest mass star which will create a core-collapse supernova 
(« 8M Q ). 



Conservation of mass 



Ein — /wTO-SNnsN 



0.4 



(13) 
(a,b) 



fSN/^osc (c) 



This equation incorporates four processes, the first of these 
being the infall of gas onto the disk, Ei n which is calculated 
from the technique described in section |2.2| 

Gas is converted to stars, hence the term — E* , and there 
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will also be mass lost from the system due to winds, a process 
parameterised by / w - Kampakoglou & Silk (2006) examine 
a range of constant values for this parameter between and 
1 but a single value / w = 0.4 is used here in the model which 
reflects their work (b), to avoid varying too many factors. 
This simple approach is also applied in model (a). 

Model (c) introduces a functional dependance for / w 
which is motivated by the following simple argument. If mo- 
mentum is conserved perpendicular to the galactic plane in 
the case of a supernova releasing energy Esn = 10 44 J and 
ejecting mass AmsN = 10M©. 



E ou t < 



V2£ SN Am SN . 80kms _1 
n S N < 



■msNiisN- 



(14) 



Relativistic corrections to this argument would be small 
due to the low speeds involved and would only serve to 
lower the effective ejected mass, further enforcing the 
inequality. If this limit does ind eed govern the process on 
an approximate leveQ equation \2.3\ can be changed to the 
equality applied in equation (13:), hence usn = 80kms -1 . 




t / Cyr 



Heating and cooling of the interstellar medium 



(15) 



^6 



EsNhsN 



Q^cEcool 



2E E 



This equation simply applies approximate conservation of 
energy to the disk of gas clouds. Energy density, Egcr 2 , is 
gained from supernova explosions, reduced by the infall of 
cold halo gas and radiated away through collisions between 
gas clouds. 

The expression for cooling rate follows (Efstathiou 
2000) and Kampakoglou & Silk (2006), with all constants 
absorbed into the parameter E coo i = IOMqpc -2 . The value 
used by Kampakoglou & Silk (2006) for the fraction of su- 
pernova energy absorbed, e c = 0.03, is also followed here. 

By making the working approximation E* « 
( 1 — / r ) msNnsN, the interesting case where heating and 
cooling rates are equal can be evaluated analytically, giv- 
ing an expression for the velocity dispersion at this point: 



-Esn 
msN 



(16) 



This is plotted in Figs. [5] & [3] and used to help understand 
the results of the full dynamic calculation in terms of this 
simple approximation. 



Figure 4. Disks formed in isolation: Star Formation His- 
tories of the systems from Figures [5] and [3] The lower panel 
simply shows the total star formation in the disk as a function of 
time under each of the three models (a), (b) and (c). The upper 
panel shows a disk-averaged value for the effective star formation 
efficiency per orbit r E*/Eg. 



3 APPLICATIONS 



Firstly, in §3.1[ the model is used to examine whether it 
might be possible to identify the predominant mechanism 
of star formation directly from a disk galaxy's final state 
(Figj3]) or indirectly from its formation history (Figs.& [4] 
[5] ). For the purposes of this initial examination, and by 
way of introducing the model as a whole, disks are formed 
within an isolated halo which does not undergo hierarchical 
growth but instead is assumed to have existed since t = 
0. The particular parameters of the halo are the same as 
were used to produce Figure [2] Since one of the ultimate 
goals of this work is to assess the importance of hierarchical 
growth on disk formation models, sections |3.2| and |3.3| have 
been included to provide a brief illustration of its effect. 
Only one of the three models (a) is used for the purposes of 
studying this additional complexity, just as consideration of 
halo growth and mergers is set aside while investigating star 
formation mechanisms. 



2 Conserving energy for this process gives the bound: 

. 2£ SN . / lOOOkms^V 

Eout < — — n SN < I m SN n SN . 

^esc \ ^csc J 

So, even if as little as 3% of this energy were actually available, 
the constraint of momentum conservation | |2.3| is still the more 
relevant for v cac < lOOOkms - 1 , which is the case for the majority 
of disks. 



3.1 An Isolated Halo 

Many key features of the three resulting systems shown in 
Figure [3] show no clear evidence of the particular underly- 
ing physics. The considerably different assumptions made, 
notably regarding gas outflow, affect the overall mass dis- 
tribution, and hence the circular velocity, by no more than 
10%. Consequently, the scalelengths of the exponential pro- 
files which best match the stellar surface densities are also 
virtually identical, since nearly all the inner part of each disk 
has been converted into stars. 

So, here are two principle observables of a disk galaxy, 
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Figure 3. Disks formed in isolation: Radial profiles of disks formed within an isolated halo with parameters as specified in the text. 
Columns correspond to the three different models of star formation described in | ]2.3| (see labels at the top of each column) . Throughout 
the Figures, dotted blue and dashed red lines represent the gaseous and stellar components respectively. Row 1 (the topmost) shows the 
gravitational stability parameter. In Row 2 the solid lines show the galaxy rotation curve. The contributions of the disk and the dark 
matter are indicated by faint dashed and dotted lines respectively. Green solid lines give the gas velocity dispersion, <r g and dotted lines 
show the value that would be expected if heating and cooling rates were equal (see text). In Row 3 the star formation efficiency is shown 
by the solid yellow curve, while the others show the metallicities of gas and stars. Row 4 shows surface density profiles. The dashed black 
line indicates the total surface density in the absence of any outflow and the solid green lines show the exponential profile which best 
matches the stellar density. Scalelengths are (a) 1.34, (b) 1.34 and (c) 1.33 kpc. 



rotation curve and surface brightness profile, which can be 
explained equally well by three competing evolutionary the- 
ories. This is an ideal platform from which to ask what fur- 
ther observations would be required to break this degener- 
acy. The required details suggested by this simple study are 
the metallicity and velocity dispersion of the older (in this 
case inner) parts of the disk. 

The reason for differing velocity dispersions is most eas- 
ily appreciated by considering energy balance in the ISM, an 
issue already alluded to in eqn. \16\ . The velocity dispersion 
which corresponds to this state is shown by a dotted line 
in Fig[3] demonstrating that the value predicted by more 
dynamic modelling (solid line) is very close to this in the 
absence of recently infalling cold gas (which occurs mostly 
in the outer parts of the disk at this final stage). The ex- 
ception to this trend is the inner parts of the disk in model 
(b), where this equilibrium has evidently not been attained, 
recent star formation still heating the region effectively. 

In models (a) and (b) here, the star formation rate, and 



hence heating rate, is assumed to be proportional to epicyclic 
frequency. The cooling rate has a similar dependance so, 
under the condition of energy balance, 

ff g ocEj*. (17) 

In the collisional model of star formation however, the heat- 
ing rate is independent of dynamical timescale and the same 
condition leads to <r g oc r c . Since the epicyclic period is 
shorter toward the centre of the disk, this explains the lower 
central velocity dispersion in Figure [3jc) . This prediction is 
interestingly at odds with observations (Kamphuis & Sancisi 
1993; Schuster et. al. 2006) which find that <r g monotonically 
decreases. 

Taking a global view of the evolution, shown in Figure[5] 
reveals that the systems cool below their equilibrium point 
after the initial burst of star formation and are kept at this 
lower level due to continuing infall of cold gas. This effect is 
most pronounced in (b). 

All three theories predict similar eventual star forma- 
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Figure 5. Disks formed in isolation: Histories of the systems shown in Figure [3] The cooling time at the virial radius is 12.5 
Gyr. Row 1 shows mass-weighted mean gravitational stability and row 2 shows the circular velocity averaged over the stellar mass. The 
velocity dispersion averaged over the mass of gas is shown by the solid line in row 3 and the dotted line again represents the value that 
would be expected from energy balance. Row 4 shows the metallicities and row 5 shows the total masses, the dashed black line again 
representing the hypothetical case of no outflow. The bottom row shows the stellar disk scale length. 



tion efficiencies. However, the predicted histories differ con- 
siderably, as can be seen in Figure [4] The clearest feature 
is that a dynamic star formation mechanism could pro- 
duce much higher initial efficiencies whilst still being con- 
sistent with the relative inefficiency which is observed at the 
present-day . 

According to model (b) , this behavior is caused by high 
instability in the embryonic disk. As the cold gas is depleted, 
the disk stabilizes and the star formation efficiency duly 
drops to its final value. The evolution of the stability param- 
eter Q, along with other quantities, can be seen in Figure 

El 

Under the assumptions of model (c), it is the frequent 
cloud collisions in the highly dense, early formed inner re- 
gions which cause the high efficiency. Subsequent loss of ef- 
ficiency as the gas density drops is exaggerated further by 
the accompanying drop in velocity dispersion (Fig[5| as su- 
pernova heating decreases along with star formation. This is 
why the third model tends more rapidly to a low efficiency. 

All these differences remain fairly subtle, suggesting 
that the final state of a galactic disk, however uniquely it 
may have been reached, is predominantly governed by gas 



supply, gravitational dynamics and stellar lifetimes, which 
are the same for all the models presented here. The effects 
of hierarchical growth, as shown in the next section, would 
swamp most signatures of the star formation process. 

Other models introduce bulge formation through disk 
instability. Van den Bosch (2001) examines the fractional 
contribution made by the disk to the circular velocity Udisk/^ 
and assumes that disks where this exceeds 0.7 at any point 
are globally unstable. The disk component of these systems 
is indeed larger than this value, but only marginally (just 5- 
10% of the mass would have to be redistributed in order to 
satisfy the criterion). Avila- Reese k, Firmani (2000) assume 
instead that stars from regions exceeding a critically stable 
surface density, E* = (j*/3.36r e G, are transfered from the 
disk to the bulge until the equality is satisfied. This would 
have a very dramatic effect on the system in this example, 
removing much of the inner part of the disk (r < 5kpc) and 
leaving a remainder with little resemblance to an exponen- 
tial profile. 

A clear set of theoretical or observational evidence is 
still needed to support one particular analytic approach to 
this issue. In lew of this, we note the theory that perturba- 
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Figure 6. Column 1— Hierarchical formation with no 
galaxy merging: See the caption to Fig. [3] for the key. The 
exponential scalelength which best fits the stellar surface density 
is 1.5kpc. Column 2 — Galaxy merging Included: . Following 
the formation of a central bulge with mass 2.9 X 10 10 Mq , the stel- 
lar disk scalelength reduces to l.Okpc. The dashed black line in 
the plot velocity vs. radius still shows just the disk's contribution 
to the rotation curve. 



Figure 7. Column 1— Hierarchical formation with no 
galaxy merging: See the caption to Fig. [5] for the key. Col- 
umn 2 — Galaxy merging Included:. The growing mass of 
the bulge component is shown here by the solid red line in the 
plot of Mass vs. time, culminating in a final Bulge/Disk ratio of 
0.45. The dotted red line still shows the total stellar mass, which 
now includes the stars in the central bulge. 



tions on scales larger than A c = 4-7T GEr c will not propagate 
even if Q > 1 (Binney & Tremaine 1987). In this particu- 
lar example, A c is never more than about a tenth of that 
region's galactic radius and is usuallly much less. 



3.2 Hierarchical Growth 

The calculation of disk formation is now embedded within 
the Galform semi-analytic framework. A statistically rep- 
resentative sample of dark matter halo merger trees was 
generated, corresponding to present-day halo masses in 
the range 10 12 to 1O 13 /i -1 M0. From these disk-dominated 
galaxy at z — has been selected which shares similar prop- 
erties with the isolated system in §3.1| These calculations 
used the standard Galform assumptions, adopting param- 
eters chosen by Baugh et al. (2005). 

A full description of the hierarchical growth of struc- 
ture, and how this affects the formation of galaxies, can be 
found in Cole et al. (2000). We will briefly discuss aspects 
of this implementation central to the present work. 



The growth of structure in the dark matter component 
is followed by constructing merger trees on a fine grid of 
timesteps. Within this merging hierarchy of halos we then 
proceed to follow galaxy formation. The first halos to form 
are assumed to have a fraction f2b/£2o of their mass in the 
form of baryons which are shock heated to the halo virial 
temperature. We then track the cooling radius within this 



atmosphere of hot baryons (see [2.1 1 to determine the rate 



at which gas cools. Two aspects of this cooling are crucial 
for the present work: 

(i) Later generations of halos will contain a smaller fraction 
of baryons in the hot atmosphere as some will have been 
locked up into galaxies at earlier times; 

(ii) The growth of halo mass is a continuous process, but 
within Galform, we impose discrete halo formation events 
which occur whenever a halo has doubled in mass since the 
previous formation event. At such events, the temperature 
of the hot gas is reset to the new virial radius and the 
cooling radius is reset to zero (beginning to grow again as 
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described in f 2.1 1 



The second point has been a traditional feature of semi- 
analytic models (originally arising for pragmatic reasons). 
It is retained in this work in order to maintain similarity to 
previously published work but its lack of physical motivation 
is acknowledged and a more realistic cooling model (such as 
that of McCarthy et al. 2007) is planned. 

In the first example here, illustrated in the first column 
of Figures [6] and [7] the effects of mergers between galaxies 
are ignored. Specifically, the effects of dynamical friction are 
not included, so sub-halos and their galaxies will continue 
to orbit within their host halo forever. As such, no elliptical 
galaxies or bulges can forn^] 

Several differences with respect to the isolated case are 
immediately apparent. Firstly, star formation now occurs in 
short bursts, even though the star formation is occurring 
purely in disks. This is due to the hierarchical build-up of 
structure — significant merger events replenishing the sup- 
ply of hot gas which is available to cool onto a forming disk. 
(This also causes a decline in the gas metallicity.) This cor- 
respondence can be seen directly by comparing the times of 
peak star formation rate with the step-like increases in disk 
mass and rotation speed. It should be noted that observa- 
tions of the Milky Way are consistent with a discontinuous 
star formation history (Rocha-Pinto et al. 2000). 

Though qualitatively similar to the rotation curves 
shown in Firmani & Avila-Reese (2000) and van den Bosch 
(2002), those predicted by this model are definitely less con- 
sistent with observations, the rotation speed being higher in 
the disk centre for reasons discussed in §2.2| 

After an initial high peak, the velocity dispersion of 
the hierarchically-formed systems tends to increase slightly 
with radius, as would be expected from the approximate 



dependance in equation (171. This is again in marginal dis- 



agreement with observations, suggesting that this evolution 
model may be too simplistic. 

The surface density profiles of the stellar disks are no 
longer completely smooth, due to the discontinuous nature 
of hierarchical growth. This problem is believed to be exac- 
erbated somewhat by the nature of the cooling model used 
in Galform. As described by Cole et al. (2000), the cooling 
time for hot halo gas in Galform is reset to zero each time 
a halo is considered to have "formed" (formation events are 
taken to occur whenever the halo's mass has doubled since 
the previous formation event). 



3.3 Galaxy Mergers 

Finally, we allow galaxy-galaxy mergers to occur and, in 
the case of major mergers, to result in the formation of 
spheroidal systems which may later become the bulges of 
disk-dominated galaxies. In the Galform model, satellite 
galaxies within a halo can lose orbital angular momentum 
via dynamical friction and eventually merge with the galaxy 
residing at the centre of that halo. The results of such merg- 



3 Note that the Baugh et al. (2005) model does not allow for 
the formation of spheroidal systems through global instabilities 
of disks. 



ers are found by applying the usual rules of Galform (see 
Cole et al. (2000) for details). 

If disks are destroyed by mergers, we use the properties 
at the half-mass radius of the disk (now determined explic- 
itly by our calculations) in determining the properties of the 
resulting elliptical galaxy. In minor mergers, disks are not 
destroyed but added together by assuming that the cumu- 
lative mass as a function of specific angular momentum is 
the same for both the combined two-disk system and the 
resulting, single disk. The new radial density profile is then 
determined as described in §2.2| 

In the presence of a bulge, eqn.[6]is modified to become 



R = 



GfMdiskt?) + MdmO') + Mbulge(-R)] 



(18) 



and is solved iteratively. Mb u i gc (-R) is the mass of the bulge 
within radius R. We assume a de Vaucouler's profiles for 
bulges, and solve for the characteristic radius using the same 
methods as Cole et al. (2000). 

This treatment neglects the orbital angular momentum 
possessed by the galaxies just prior to merging. We are not 
aware of any detailed (e.g. N-body) experiments which have 
assessed how this angular momentum is redistributed into 
the resulting galaxy (or if, instead, it is mostly lost from the 
system), although a merger-driven disk formation scenario 
has been proposed by Robertson et al. (2006). 

To quantify the importance of this orbital angular mo- 
mentum, we performed a simple test. Specifically, prior to 
merging the satellite galaxy onto the primary galaxy we in- 
cresed the specific angular momentum of each annulus of 
material in the satellite disk by an amount equal to its or- 
bital specific angular momentum. We found essentially no 
change in the resulting profiles of disks at z — and there- 
fore conclude that a correct treatment of this orbital angular 
momentum would likewise cause only minor changes in the 
predicted galaxy profiles. 

The resulting galaxy disks shown in the second col- 
umn of Figures [6] and [jj are very similar to the mergerless 
case, despite the fact that approximately 30% of the final 
galaxy mass is in a bulge component. One significant and 
expected difference is in the rotation curve which rises even 
more rapidly towards the centre due to the additional cen- 
tral mass. 



4 CONCLUSIONS 

We have described methods for expanding the range of phys- 
ical properties predicted by semi-analytic models of galaxy 
formation to include the surface density profiles of stars 
and gas in galaxy disks. These profiles are computed self- 
consistently from the assumed distribution of angular mo- 
mentum in the hot halo gas which fuels galaxy formation. 

Using the additional information provided by these re- 
solved disks we are able to explore more detailed models 
of star formation than has previously been possible within 
semi-analytic models. One key difference is the provision for 
dynamic evolution of the system, free from any prior en- 
forcement of fixed stability (Firmani & Avila-Reese 2000; 
Efstathiou 2000) or energy balance (van den Bosch 2001; 
Kampakoglou & Silk 2006). 

These methods have been implemented within the Gal- 
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FORM model in order to study how disks form and grow in 
a hierarchical Universe, exploring three cases of increasing 
complexity and realism. 

The main conclusion from this initial, brief study is that 
the primary factor governing the evolution of these systems 
is the availability of gas for star formation, the introduction 
of hierarchical growth dramatically altering the form of the 
resulting star formation histories. This emphasizes the fact 
that any realistic model of disk formation must be embedded 
within the appropriate hierarchically growing distribution of 
dark matter halos. 

By comparison, the structure and internal properties of 
the final disk systems are relatively insensitive to the as- 
sumptions concerning star formation. The chemical evolu- 
tion and star formation histories exhibit greater distinction, 
so this conclusion can not necessarily be extended to galax- 
ies of all masses at all redshifts — further work is needed to 
explore these regimes. 

There is no doubt that a prescription for star formation 
could be created which did produce contrasting results, but 
the fact that three different, physically motivated star for- 
mation models produce quite similar present-day systems 
suggests that semi-analytic models (and hydrodynamical 
simulations which frequently use similar rules for the sub- 
grid star formation physics) can produce robust predictions 
for galaxy properties even though the precise rules governing 
star formation are not known. 

Another key result is that current estimates of the dis- 
tribution of angular momentum in halo gas from up-to-date 
numerical simulations (Sharma & Steinmetz 2005) corre- 
spond to an excessive quantity of low angular momentum 
material, consistently producing disks whose rotation curves 
rise too sharply in the centre, a problem which has been well 
noted (Bullock et al. 2001; van den Bosch, Burkert & Swa- 
ters 2001; van den Bosch 2002; Mailer & Dekel 2002). Even 
more accurate characterisations of the distribution of halo 
gas angular momentum, and further studies of how this is 
conserved (or otherwise) during cooling and star formation 
are keenly awaited. 

While this work explores the evolution of a single sys- 
tem, it is just one of a large, statistically representative sam- 
ple of galaxies forming in halos of a wide range of masses. 
The full potential of this modeling technique will be ex- 
ploited in forthcoming work, its predictions being set against 
both cosmological data-sets and collected galactic profiles 
(of stellar and gas mass, rotation speed, age and metallic- 
ity). This will provide an important assessment of the com- 
patibility between ACDM cosmology, as represented by its 
predictions for large scale structure formation, our under- 
standing of star formation and the observed Universe. 
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APPENDIX A: NUMERICAL ROBUSTNESS 

As described in §2.2| processes internal to each disk are fol- 
lowed for each of 200 concentric annuli and external pro- 
cesses, e.g. galaxy mergers, are followed on a coarser set of 
20 "zones" (also concentric annular regions, but a different 
term is used to avoid confusion). 

To establish the numerical robustness of these approx- 
imations, we have repeated our calculations using different 
numbers of annulli and zones. Figure |A1| shows the effect 
these alterations had on the eventual profile of the typical 
galaxy examined in this paper and Figure |A2| shows the 
fractional changes in each parameter. 

It is representative of the finding, for all of the modelled 
galaxies, that increasing the number of zones above 20 made 
no more than a 5% difference in any eventual parameter 
value within the inner 90% of the disk's mass. In general, 
the difference was much less than this. Our calculations are 
therefore relatively unaffected by these numerical choices. 
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Figure Al. Left panel: Radial profiles of the galaxy examined in 
Sounder model (a), as generated using different numbers of radial 
divisions used in calculations of mergers ("zones"). Lines have 
the same meaning as in Fig. [3] The disk mass and gas fractions 
given in the panels correspond to the number of Zones with the 
corresponding typeface. The vertical dotted line encloses 90% of 
the disk's mass. Right panel: Fractional differences between each 
pair of lines. 



Figure A2. As for Figure |A1| but now varying the number of 
radial divisions used in calculations concerning the internal evo- 
lution of each galaxy ("annuli"). The vertical dotted line again 
encloses 90% of the disk's mass. 



